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We present a detailed experimental procedure for preparing relativistic vortices, governed by the 
nonlinear Dirac equation, in a two-dimensional Bose-Einstein condensate (BEC) in a honeycomb 
optical lattice. Our setup contains Dirac points, in direct analogy to graphene. We determine a 
range of practical values for all relevant physical parameters needed to realize relativistic vortices 
m BEC of ^'^Rb atoms. Seven distinct vortex types, including Anderson- Toulouse and Mermin-Ho 
skyrmion textures and half-quantum vortices, are obtained, and their discrete spectra and stability 
properties are calculated in a weak harmonic trap. We find that most vortices are stable with a 
lifetime between 1 and 10 seconds. 



Multi-component Bose-Einstein condensates (BECs) 
present an ideal setting for studying complex vor- 
tex structures [1]. Such vortices allow for topologi- 
cally intriguing configurations ranging from skyrmions 
to knots [2|-y. The usual method for adding a spinor 
structure to a BEC relies on hyperfine degrees of free- 
dom or different atomic species. Instead, we use the band 
structure and linear dispersion relation around the Dirac 
points at the Brillouin zone edge of a honeycomb optical 
lattice to realize a four-component Dirac spinor, in direct 
analogy to graphene This gives us both pseudospin 
as well as a relativistic structure. To accomplish this, we 
propose starting with a BEC of weakly interacting alkali 
metal atoms in the lowest Bloch state of a 2D honeycomb 
optical lattice, then using Bragg scattering to populate 
Bloch states at the two inequivalent Dirac points, fol- 
lowed by the application of a Laguerre- Gaussian laser 
beam to deliver a net angular momentum to the BEC 
which excites a plethora of vortex structures. The vor- 
tices we obtain are solutions of the nonlinear Dirac equa- 
tion (NLDE), whose stability is determined by the rel- 
ativistic linear stability equations (RLSE) [HI, Our 
work on the NLDE+RLSE system opens up the field of 
relativistic simulations in BECs at velocities ten orders 
of magnitude slower than the speed of light. 

In this letter we combine the study of Dirac points with 
superfluid vortices, an environment reminiscent of par- 
ticle physics models where relativistic vortices are com- 
monplace isl^llll- Stability of a BEC at the Dirac points 
presents a challenge, since Bloch states there have fi- 
nite crystal momentum and nonzero energy. We handle 
this problem by introducing an intermediate asymmetry 
between the A and B-sublattice potential depths which 
opens up a mass gap. Using a gap enables us to con- 
struct initial and final Bloch states, i^A.o and V^a,k (with 
Dirac point momentum K), as superpositions of the two 
degenerate states at the Dirac point with velocities ci 
and — Q, respectively. This produces a state with group 



velocity equal to zero, relative to the lattice. Stationar- 
ity of the BEC with respect to the lattice and the lab 
frame is significant experimentally, since the BEC can 
remain confined in an external trapping potential indef- 
initely and does not suffer from dynamical instabilities 
associated with relative motion between the BEC and 
lattice. The end result is a metastable state in which 
thermal losses can be managed by maintaining the sys- 
tem at very low temperatures. For realistic experimental 
parameters our relativistic vortices are stable for up to 
10 seconds, as long or longer than the lifetime of typical 
BECs. 

Relativistic vortices are realized in the emergent non- 
linear Dirac background, in the long wavelength limit 
of a 2D honeycomb lattice. The usual BEC parame- 
ters in 3D are renormalized, once for the dimensional 
reduction [l8| , and again after integrating over the lattice 
Wannier functions and going to long wavelengths. Conse- 
quently, NLDE physics is only experimentally realizable 
in practice when several energy and length constraints are 
satisfied. We list these constraints in Table |T] along with 
their mathematical definitions. For our calculations, we 
use the semiclassical estimate [l9| of the hopping param- 
eter th = 1.861 (Fo/^i?)'^^^i?exp(-1.582yV^)- It 
is helpful to consolidate the constraint inequalities to ar- 
rive at expressions relating the temperature T and length 
scales of the system, a^, a, L^, and R^: 



1 < 



3/2 



< 



2^v/2^V^(d^a,)V^ 
3V3a2 [l + 7ra/(4y2i?x)] ' 



(1) 
(2) 



where d is the average inter-particle distance defined in 
terms of the particle density d = n~^/^. All other quanti- 
ties are defined in Table HI The temperature T in Eq. ([2]) 
depends indirectly on the ratio Vq/Er through n. We 
can get an idea of how the particle density affects T by 
evaluating the inequalities for different values of n while 
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Parameter 


Symbol/Definition 


Constraint 


Value 


Range 




(a) Temperature 


T 




^ hwz 


8nK 


- 100 pK <T < 


- 80 nK 


(b) Chemical potential 


fi 




^ hwz 


2.36 nK 


< 4.10 nK 




(c) Transverse oscillator length 


T 


/fr /7\/f, , \l/2 


< R± 


3.0 /im 


< 2.10 /im 




(d) Healing length 


^ = 


- l/\/87rnas 


^ Lz 


2.14 /im 


< 1.66 /im 




(e) Effective speed of light 


ci -- 




< Cs,2D 


2.69 X 10"^cm/s 


< 5.40 X 10~^cm/s 


(f) Dirac nonlinearity 


U : 




< th, M 


0.393 nK 


< 2.36 nK 




(g) Quasi-particle momentum 


k = 


--p/h 


< V8/a 


6.27 X 10^ cm"^ 


6.27 X lO^cm"^ 


< < 5.66 


(h) Dirac healing length 




rac = thaV^/2U 


> a, < 7?^ 


5.25 fim 


0.50 /im < Coirac 


< 50.0 /im 


(i) Lattice depth 


Vo 




> Er 


2.59 fiK 


0.79 fiK<Vo <l 


S.95 fiK 



TABLE L Physical parameters and constraints for the NLDE typical for a BEC of Rb atoms. (a,b) Relative energies for 
the 3D to 2D dimensional reduction, with the vertical trap oscillator energy hwz- (c,d,h) Relative lengths for the 3D to 2D 
dimensional reduction, (e) Landau criterion imposed to avoid dynamical instabilities, where the 2D speed of sound in the 
continuum Cs,2D = \/3gn/2M = 2.97 x 10~^cm/s. (f) The weakly interacting and superfluid (not Mott insulating) regime, 
(g) The linear Dirac cone approximation which requires that quasi-particle momenta hk remain small compared to the Dirac 
point momentum hK. (h) Long- wavelength limit, which sets the scale for the 2D Dirac healing length, (i) The lowest-band 
and tight-binding approximation. For the values in the table, we use the ratio of lattice depth to recoil energy Vq/Er = 16, 
lattice constant a = 2A/3 = 0.55 /im, and planar trap radius R± = 100 a, average particle density n = 1.5 x 10^^ m~^, hopping 
energy th = 4.31 nK, and atomic mass of ^^Rb. 



fixing Vq/Eji = 16. For example, n = 10^^ m ^ gives 
26.259 /im < < 86.934 /im and T < 8.17 x IQ-^nK, 
whereas for n = 10^^ we find 0.263 /im < < 
0.187 /im and T < 162 nK. From this we see that a prac- 
tical value for T requires that densities be considerably 
larger than 10^^ m~^, a consequence of the additional 
constraints in Eqs. ([I])-(|2]). 

The honeycomb lattice is composed of two degenerate 
hexagonal lattices, A and B, which leads to the band 
crossing and the Dirac cones. The lattice potential is 
straightforward to implement experimentally [2o|, Elj us- 
ing light tuned either to the blue or to the red of an 
atomic resonance. For the red-detuned lattice, all three 
laser fields are polarized parallel to the plane of propaga- 
tion so that the polarization of the net field is spatially 
dependent. This polarization gradient produces a spin- 
dependent lattice potential which is key to our method, 
allowing for differentiation between different hyperfine 
ground states. Figure [TJa) shows the optical potential 
produced for ^^Rb atoms in different hyperfine states 
when the lattice is formed from = 830 nm light [ioj . 

Preparation of the BEC at a Dirac point is accom- 
plished by first condensing into the lattice via evaporative 
cooling, then inducing Bragg scattering between crystal 
momenta and K using auxiliary laser fields. We start 
with the BEC in the lowest-energy Bloch state and hy- 
perfine state rriF ^ with the spin-dependent potential 
turned on so that a mass gap is in place. The gap is 
defined as 2 \ms\cf^ which measures the difference in po- 
tential depths between the A and B-sublattices as shown 
in[2fa). Note that rrisCi is the corresponding mass term in 
the NLDE, with ci the effective speed of light (see Table 
I). The lattice depth is then increased adiabatically (ioj . 
Only the sublattice with the lowest energy is occupied at 



this point, assumed here to be the A-sublattice. Bragg 
scattering to a Bloch state at a Dirac point can be ac- 
complished by applying two laser fields with wavevec- 
tors kbi and kb2, which satisfy kbi — kb2 = K and 
have frequencies cjbi and a;b2- These frequencies satisfy 
^bi — ^b2 = Ao; = [EaCK) — EA{0)]/h. Here, the func- 
tion EaO^) gives the dispersion relation for the lower (A) 
band. Application of this potential results in Rabi oscilla- 
tion between the initial state ^^^,0 (ground state), and the 
final state '0a,k (Dirac point) with a Rabi frequency I^b, 
where i/ja, V^b are wavefunctions; is the BEC 

density distribution on the A (B) sublattice. Figure [T]^b) 
shows numerical calculations for IHQbI as a function of 
the depth of the honeycomb lattice Vsc in units of the 
depth of the Bragg scattering lattice Vb - We should keep 
in mind that Vb <C Er. The entire population of atoms 
in state V^a,o can be transferred to ipA,K by applying the 
Bragg scattering potential for a time r-j^ = tt/Q provided 
that the amplitude of the Bragg potential is chosen such 
that H/tt^ is significantly smaller than the energy split- 
ting between the upper and lower bands, i.e., 2\ms\cf. 
Note that because of the gap, all the atoms are presently 
in a state with zero group velocity. 

With the BEC at the Dirac point but only in the 
A-sublattice sites, we populate the B-sublattice by first 
transferring atoms to a hyperfine state that does not ex- 
perience the vector light shift and therefore has no mass 
gap. We can then populate the B-sublattice Bloch state 
by applying a periodic perturbation which modulates the 
amplitude of one of the lattice laser fields. This provides 
an anisotropy in the tunneling matrix elements, which 
results in a net transfer of atoms to the B-sublattice, as 
depicted in both panels of Fig. [2l Finally, the transfer 
of atoms to the non-equivalent Dirac point can be 



3 




A B 



0.50 
0.45 
> 0.40 
^ 0.35 
0.30 



1 1 1 1 




-- mF = ±1 - 


1 1 


1 1 




FIG. 1: (color online) Bragg scattering in a spin- dependent 
honeycomb lattice, (a) Honeycomb potential for ^^Rb atoms 
in state |F, m^) = |2, 1) for the case when the wavelength 
of the lattice light Xl = 830 nm. (b) Rabi frequencies for 
the mp = i2 (solid blue) and dil (dashed red) states, (c) 
Rabi frequency for transitions between non- equivalent Dirac 
points for cases where the sub-lattice index remains the same 
(solid blue) or changes (dashed red) as functions of the depth 
of the scalar part Vsc of the optical lattice potential. (Inset) 
Time dependence of the sublattice populations at the Dirac 
points K and for an optical lattice depth of 14c = 4:Er. 



accomplished again by Bragg scattering from a lattice 
formed using auxiliary laser fields 



22|. The Rabi fre- 



quency and resulting time-dependent populations for the 
K and points are depicted in Fig. [TJc). 

We analytically and numerically obtain seven phys- 
ically distinct NLDE vortex types as follows (a table 
detailing the functional form of each vortex type is in- 
cluded in supplemental materials, as well as a brief re- 
view of the NLDE and RLSE). (i) The vortex/soliton is 
a bright soliton or density peak in the center in the first 
component with a vortex of phase winding 27r around 
the outside in the second, (ii) The ring- vortex/soliton is 
also a bright soliton in the first component, but the vor- 
tex component is a ring peaked near the healing length 
r = ^Dirac- (iii) The Anderson- Toulouse skyrmion has 
the same core structure as the vortex/soliton, but the 
spinor components are continuously interchanged as the 




FIG. 2: (color online) Coherent transfer between sublattices A 
and B. (a) Three step process of exciting atoms from the A- 
sublattice with hyperfine state |1,0) (no sublattice asymme- 
try) to the hyperfine state |2, 1) via the rf/mw transition mwi, 
then from the A-sublattice to the B-sublattice via the pertur- 
bation Hm, and finally back to the |1,0) hyperfine state via 
the mw2 transition (mw = microwave), (b) The same process 
showing the associated rf/mw frequencies and the detunings 
Amwi and Amw2. 



distance from the core increases, while staying within the 
bounds < IV^aI, \'^b\ < 1 and conserving total density 
IV^Ap + IV^Bp = 1- (iv) The Mermin-Ho skyrmion again 
has similar behavior near the core but the soliton (vor- 
tex) amplitude decreases (increases) monotonically away 
from the core within the bounds cos(7r/4) < ipA < 1 
and < ipB < cos(7r/4). (v) The half-quantum vor- 
tex or semion is characterized by a phase discontinuity 
such that far from the core the amplitudes have the form 
i/ja oc cos(6>/2) and i/jb oc sin(6>/2); the additional tt phase 
is accounted for by a rotation between the Dirac spinor 
components. So far, all of these solutions have one unit 
of angular momentum, £ = 1 with ^ G N, either a phase 
winding of 27r in one component or a winding of tt in 
each component. Additionally, for arbitrary phase wind- 
ing (^ > 1) (vi) ring- vortices and (vii) topological vortices 
exist with i — 1 {£) units of winding in the first (sec- 
ond) spinor component, but differ in their asymptotic 
form. Component amplitudes for the ring-vortex peak 
at around one healing length from the core and quickly 
decay for large r. On the other hand, topological vor- 
tices retain non-zero density far from the core. Several 
representative vortices are plotted in Fig. [3l All of the 
vortices here can be created using straightforward vari- 
ations of the transition sequence depicted in Fig. [2l In 
the supplemental materials, we explain how vortices are 
realized in practical experiments. 

For most of vortex types (i)-(vii), we find lifetimes 
r to be long compared to the lifetime of the BEG 
itself. Using the RLSE (as detailed in supplemen- 
tal material), we obtain we obtain the following val- 
ues for r: 9.13 s, 10.43 s, 11.51s, 1.57 x 10"^ s, 1.57 x 
10"'^ s, 1.25s, 1.29xl0~^s; for the vortex/soliton, ring- 
vortex/soliton, Anderson- Toulouse, Mermin-Ho, half- 
quantum, i = 2 ring-vortex, and i = 2 topological vor- 
tex, respectively. To complete our description, we have 
computed the discrete spectra pertaining to the case of 
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(b) arg(^)/7r 




FIG. 3: (color online) Plots of relativistic vortices. Density 
and phase of (a,b) i — 2 ring-vortex, (c,d) B-sublattice of 
Mermin-Ho skyrmion, (e,f) ring-vortex/soliton, (g,li) half- 
quantum vortex, or semion. All these vortices and more 
can be made by variations on the experimental techniques 

of Figs, na 



planar confinement by a weak harmonic trap. These are 
plotted in Fig. |4j Note that ring- vortices are minimally 
affected by the presence of a weak trap, since they are 
highly localized objects and lie very near the center of 
the trap. 

In conclusion, we have described in detail a method 
for constructing a stable BEG at the Dirac points of a 
honeycomb optical lattice. Our system allows for rela- 
tivistic vortex excitations in a macroscopic Dirac spinor 
wavefunction, providing a means of studying high energy 
field theoretic vortices in a condensed matter setting. We 
have completely specified the required physical parame- 
ters, lifetimes, and spectra for harmonically bound vor- 
tices as a prescription guide for the experimentalist. Vari- 
ations on the NLDE have tremendous potential for a host 
of relativistic simulations in BEGs. Interesting examples 
include Soler models [23| and the extended Gross-Neveu 
model [2J . Our work puts such efforts on a solid exper- 
imental footing. 
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FIG. 4: (color online) Spectra for relativistic vortices con- 
fined in a harmonic potential. (a) Vortex/soliton (black curve), 
Anderson- Toulouse skyrmion (red), Mermin-Ho skyrmion 
(blue), and half-quantum vortex (green), (b) Topological vor- 
tices for = 2, 3, 4 (black, red, blue), (c) Radial ground state 
and first two excited states of the vortex without skyrmion 
symmetry (black, red, blue). In each figure, the renormal- 
ized chemical potential is plotted as a function of the nor- 
malization (both quantities are covered in more detail in the 
supplemental material). There are two regimes characterized 
by power law: jl oc J\f^ . The weakly interacting free-particle 
regime occurs for small A/", whereas the strongly interacting 
vortex regime is in the region of large A/"). Note that the 
vertical and horizontal axes labels are dimensionless. 
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SUPPLEMENTAL MATERIALS 

NLDE Solutions - The NLDE treats the entire Dirac 
four-spinor. In its simplest realization without mass 
gaps and in tight binding the upper two components, 
called a Weyl spinor^ are decoupled from the lower two, 
and can be written ^ = {iPa^'^b)^ • We obtain vor- 
tex solutions by expressing the spinor components in the 



form: iljA{r,0,t) = ±i fA{r)e 



iii-l)6 



lie , 



-ijjit/h 



iPB{r,e,t) = 



and writing the NLDE in plane-polar 



fB{r)e 
coordinates: 

- hci (dr + l^Bir) + U \fA{rffA{r) = Mr) (3) 

^i(dr + ^yAir) + U\fBir)\^fB{r) = /x/bW, (4) 

where i is the integer phase winding and the other pa- 
rameters are defined in Table 1 of the main article. For 
the case /i = 0, Eqs. ©-(jl]) give closed form expressions 
for the radial amplitudes /a and fs - These are the ring- 
vortex/soliton = 1) and general ring- vortex (^ > 1) 
solutions. For the case /i 7^ 0, closed form solutions exist 
in some cases while others are obtained using a numer- 
ical shooting method. Solutions of the NLDE are listed 
in table HIl 

Experimental Realization of Vortices - Relativistic vor- 
tices can be excited by starting with all the atoms 
in the A-sublattice at the Dirac point, then applying 
co-propagating Gaussian and Laguerre- Gaussian laser 
beams, where the Laguerre- Gaussian beam carries a sin- 
gle unit of orbital angular momentum. The transfer of 
angular momentum to the atoms occurs via the stimu- 
lated Raman transition [25| . The spatial variation of the 



beam results in mainly the B-sublattice being populated 
(the vortex) throughout most of the 2D lattice, except 
within a small disk which becomes the core of the vor- 
tex. On the other hand, the A-sublattice is left depleted 
everywhere (the soliton) except near the core of the vor- 
tex. This describes excitation of the vortex/soliton or 
Anderson- Toulouse skyrmion [6[. The Mermin-Ho vor- 
tex can be obtained by the same process, but by only 
partially transferring atoms to the B-sublattice. The sub- 
lattice amplitudes far from the vortex core are tuned to 
satisfy IV^bP = iV^Ap < 1, where is the density 

of the EEC in the first (second) four-spinor component in 
the NLDE, and v^(^) = {h / M)V (j) a{b) is the associated 
relativistic fluid velocity, with (j)A{B) = ^^^{i^A{B)) the 
phase. The half-quantum vortex or semion can be excited 
by using a fractional optical vortex beam in order to pro- 
General 



26|,l27| 



vide the required angular phase jump 
topological vortices have phase winding ^ > 1, non-zero 
chemical potential /i, and satisfy IV^^I, IV^bI 7^ far from 
the center of the trap. These can be created by first ro- 
tating the lattice to excite the desired £—1 state in the A- 
sublattice, then following up with the Laguerre-Gaussian 
transition to produce the correct i versus i — 1 wind- 
ing differential in the spinor components. We note that 
there is an alternative approach to the rotating lattice 
method. General topological vortex excitations may be 
induced by subsequent applications of a two-photon tran- 
sition with co-propagating Laguerre-Gaussian/Gaussian 
beams which transfer the condensate between m = 
states (i.e. from F = l,m = 0toF = 2, m = Oor vice 
versa). Each two-photon transition changes the orbital 
angular momentum of both the A and B-sublattices by 
the orbital angular momentum carried by the Laguerre- 
Gaussian beams, while maintaining the desired winding 
differential between the A and B-sublattices. Finally, 
ring- vortices, characterized by /i = and IV^aUV^bI = 
far from the center of the trap, can be obtained from the 
other vortices by inducing depletion of the BEG from the 
outer edge of the trap towards the core. More details re- 
garding solutions of the NLDE may be found in Ref. Q . 

RLSE Solutions - The RLSE form a relativistic gen- 
eralization of the Bogoliubov-de Gennes equations anal- 
ogous to the relationship between the NLDE and nonlin- 
ear Schrodinger equation. Thus, in the RLSE the quasi- 
particle amplitudes u and v are each vector in form, to 
match the four-spinor (two-spinor at one Dirac point) 
they perturb from. The RLSE can be expressed in 2 x 2 
matrix-vector form: 



#*Vk - /7§Uk = -^kVk , 



(5) 
(6) 



where & and ^ are 2x2 matrices which contain the 
first-order derivatives {dx^idy) and the background BEG 
components V^a, and ^k is the 2x2 eigenvalue ma- 
trix. Note that U is the particle interaction. When bro- 
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Vortex type 



Winding 



Analytic form of ^(r) 



Topology 



Vortex/soliton i ■ 

Ring-vortex/soliton i - 

Anderson- Toulouse skyrmion £ - 

Mermin-Ho skyrmion £ - 

Half-quantum vortex £ - 



1 



Ring- vortex ^ = 2,3,4,. 

General topological vortex £ = 2,3,4, . 



-iO (r/rp) 



Vl+(r/ro)2 
9 (r/rp) 



Vl+(r/ro)2 



Vl+(r/ro)4' ^l+(r/ro)4_ 

ie~'^^cos(p{r/ro), smip{r /ro)] 
'^cosv?(r/ro), sin(^(r/ro)] 

T 



\ze 



[icosO/2, sin6>/2]' 



Numerical shooting method 



I^aMI = 1 

non-topological 

(p{oo) = 



: 7r/4 
= 1 



non-topological 

|^a(oo)| = 1 



TABLE II: Vortex solutions of the NLDE. Solutions are described by their phase winding, closed-form expression, and topological 
properties. Solutions which retain non-zero density far from the core have an associated conserved topological charge, and we 
state their asymptotic form. Note that ro is the length scale associated with the chemical potential or the interaction strength 
depending on the particular solution. 



ken down, Eqs. ([5])-(|6]) form a 4 x 4 eigenvalue problem in 
the quasi-particle amplitudes U]^^a{b) and V]^^a{b) (with 
momentum k) associated with particle and hole excita- 
tions of the A(B)-sublattices at a Dirac point. Vortices 
possess cylindrical symmetry so we express Eqs. (|5j)-(|6j) 
in plane-polar coordinates, factor the quasi-particle am- 
plitudes into radial and angular parts, then substitute 
in the particular solution for V^a(b)- We then obtain 
a set of first-order coupled ODE's in the radial coordi- 
nate to be solved consistently for the functions ua{b)(t)i 
^A{B)(j) and the associated eigenvalues. We discretize 
the derivatives and functions using a forward-backward 
average finite-difference scheme, then solve the resulting 
discrete matrix eigenvalue problem using a standard nu- 
merical diagonalization method. 

Vortex Lifetimes - To compute vortex lifetimes, we 
solve the RLSE to obtain the quasi-particle spatial func- 
tions and eigenvalues. In general, for vortex solutions of 
the NLDE certain eigenvalues and eigenmodes key to un- 
derstanding the physical motion correspond to Nambu- 
Goldstone modes, i.e., anomalous with a small imaginary 
component [l^. When thermal losses are small, it is the 
imaginary part of the linear eigenvalues which depletes 
the BEG. We define the vortex lifetime by computing the 
time for depletion to reach a significant fraction of the 
total fixed number of atoms in the system, and consider 
only depletion coming from the mode with the largest 



imaginary term in its eigenvalue. The lifetime is then 
given by r = [h/l-m{E)]hi{R^/ 1)^ expressed in terms 
of the largest linear eigenvalue E and the planar radius 
of the BEC i?^, in units of the lattice constant a (see 
Table 1 in the main article). Note also that the spatial 
integral / here is specific to each vortex type. For the 
experimental parameters of Table 1, we find the longest 
lived solutions to be the vortex/soliton and Anderson- 
Toulouse vortex with r = 11.51 s, compared to the typi- 
cal lifetime of a ^^Rb condensate in an optical lattice of 
less than a second [i^. 

Chemical Potential Spectra - In order to have a clear 
comparative prediction for energies involved in creating 
our vortices, we solve the NLDE using a numerical shoot- 
ing method in the presence of a weak harmonic trap of 
frequency uo^ = 27r x 0.0387 Hz along the direction of the 
lattice. This is the frequency associated with a planar 
BEC radius equal to 100 times the lattice constant. In 
this case vortices come in radially quantized states. For 
simplicity, we focus mainly on the lowest radial excita- 
tion. Using a generalization of the method in ^o|, we 
have obtained the dimensionless (renormalized) chemical 
potential /i = fi/huj_\_ as a function of the normalization 
JV = \^hw±NU/3t'f^ for each vortex type, as shown in 
Fig. 4 in the main article. Here, N is the number of 
atoms in the system with the other quantities defined in 
Table 1. 



